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Abstract 

We propose a method of calculation of the power spectrum of cos- 
mological perturbations by means of a direct numerical integration of 
hydrodynamic equations in the Fourier space for a random ensemble of 
initial conditions with subsequent averaging procedure. This method 
can be an alternative to the cosmological N-body simulations. We 
test realizability of this method in case of one-dimensional motion of 
gravitating matter pressureless shells. In order to test the numeri- 
cal simulations, we found an analytical solution which describes one- 
dimensional collapse of plane shells. The results are used to study a 
nonlinear interaction of different Fourier modes. 
Key words: large scale structure, cosmological perturbations, hydro- 
dynamics. 

1 Introduction 

Theoretical investigation of the cosmological inhomogeneity growth presents 
one of the most serious challenges to computational astrophysics. A num- 
ber of problems arises on non-linear scales when the density contrast cannot 
be assumed small. These include simulations of the galaxy formation pro- 
cess, working out predictions concerning the galactic environment (number of 
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dwarf satellite galaxies) and structure of their central regions (the cusp-core 
problem) either with cold dark matter (DM) or within warm DM models 
[H |2j [3J IU [5j [6] . Interesting possibility to obtain bounds on masses of DM 
particles stems from observations of Ly-a forest (see, e.g., |?j) and refer- 
ences therein), which requires accurate calculations of the power spectrum 
of cosmological inhomogeneity on kiloparsec scales. 

Most developed computational techniques to study the cosmological struc- 
ture formation involve N-body simulations combined with the smoothed par- 
ticle hydrodynamics [8j HI [ID] . Currently performed simulations involve up 
to 10 9 particles (see, e.g., [HE]). These methods should be tested in inde- 
pendent simulations. 

On the other hand, some analytical and semi-analytical schemes were 
proposed to study the matter power spectra on small scales [TU [TJl [13]. 
These typically involve perturbative schemes after transition to the Fourier- 
transformed hydrodynamical variables. There are techniques dealing with 
correlation functions in the Fourier space p3l [15] . These authors use some 
additional apriori suggestions in order to close the infinite chain of correlation 
functions. The validity of these suggestions is not evident. A comparison of 
different approaches can be found in [16]. 

In this paper we propose an alternative method, which uses a direct inte- 
gration of hydrodynamical equations in the Fourier space. In order to esti- 
mate workability of the method, as a first step we consider a one-dimensional 
problem of hydrodynamical evolution for a pressureless gravitating matter, 
i.e. one-dimensional density shells. The integration is performed for each re- 
alization from a random ensemble of initial data with subsequent averaging 
procedure. In order to test the numerical simulations, we found an analytical 
solution which describes one-dimensional collapse of the plane shells (Section 
[2]). This solution is used in case of periodic initial data corresponding to a 
symmetric motion. In section [3] we write down the equations for the Fourier 
coefficients and present some results for the power spectrum obtained after 
the statistical averaging. 

2 Implicit analytical solution 

In order to test numerical simulations, it is useful to have an exact solution. 
In this section we obtain such a solution of one-dimensional problem in the 
Lagrange variables. The 1-dimensional version of hydrodynamical equations 
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(continuity, Euler and Poisson equations) in case of a pressureless gravitating 
fluid is: 

— + 7— -— (2) 

dt dx dx ' 

V is the velocity, p is the mass density, $ is the gravitational potential, G 
is the gravitational constant. In order to pass to the Lagrangian description 
we introduce the stream lines X(£,t): 

— = V(X,t), X&0) = £ . (4) 

Equation (j2J) yields 

- W = E { X, t ), *<*) = -_ (5) 

Let for t = 0: p(x,0) = p (x), V(x,0) = V (x). Using (□), (EJ we see that 
£), t) — E(X(^ 2 ,t),t) = const is constant along the stream lines. 
This is a consequence of mass conservation in plane layer between X(£i,t) 
and X(£ 2 ,t). Then we obtain because of (jSJ) 

6 

X{Z u t)-X{Z 2 ,t)=Z 1 -& + t [Vo{Ci)-V (^)}-27rGt 2 J dxp (x). (6) 

Formula allows to get a general solution of the one- dimensional problem 
in Lagrangian coordinates. Further for simplicity we deal with a symmetric 
motion (with respect to the coordinate origin £ = 0) so that X(0,t) = 0, 
Vb(0) = 0. In this case 

£ 

X(^t)=£ + tV (O-2nGt 2 J dxp (x) (7) 
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The solution is valid during a limited time until 
dX 



^ = l + tV^)-2nGt 2 p (O ^0 (8) 
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Transition to Euler variables is 



p(x,t) = p (0 

where £ = £(x, t) is defined implicitly as the solution of the equation: 



V{x,t) = VoiO-AnGt j dxpo(x) (9) 

o 



x = X{£,t) 



(10) 



In case of homogeneous initial conditions p(x, 0) = pho{x) = const, V(x, 0) = 
Vho(x) = H x, H is the "Hubble constant", we get: 



x = £ (1 + H t - 27rG Ph0 t 2 ) , Ph (x, t) = - 



PhO 



H t - 2nGp ho t 2 



V h (x,t) = x 



H - AnG t p h Q 
1 + H t- 27iG Ph0 t 2 



(11) 



or 



p h (x,t) = -ggr, V h {x,t)=x^§-, R(t) = l + H t-27tGp h ot 2 . (12) 



R(t) 



R(ty 



To consider deviations from the homogeneous background, we impose 
periodic boundary conditions such that the initial conditions (t = 0) are as 
follows: 



V (0 = Ho 



oo o 



n=l 



Po(0 = PhO 



l + ^&°cos(& n £) 



n=l 



where L is the periodic "box" size, k n = 2imjL. 

Substitution of the initial conditions into (JTH9J) yields 



oo o 

az 



x(t, t) = R(t)t + tH J2r sin ^) - z (* H ^ E v sin ^) ' 



n=l 



n=l 



p(x,t) 



1+ E6°cos(fc„0 

n=l 



R(t)+tH £ aOcos^) - ! (^o) 2 E &n cos(fc n O 



n=l 



n=l 
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The solution shows a singular growth (collapse of one- dimensional gravi- 
tating layers) as the condition flS]) is violated. 

At the end of this Section we note that in spite of the exact form of the 
solution given by ([319]) we cannot avoid a numerical work when we pass to 
the Euler variables and then to the Fourier representation. So we refer to 
the method of this Section as " semi-analytical" . 



3 Numerical simulations 

Furthermore for the homogenous background we denote H(t) = dR/dt, r is 
the "conformal time": dt = R{r)dr, and y = x/R(r) is the comoving spatial 
coordinate. After some calculation on account of f[T2"j) we get 



1 + 3H t 
2H "" [3(2 - Hot) 



-In 



„ , , . dR 1 dR „ , , „ . 



Hereafter 5 is the density contrast, 9 = dv/dy, v is the peculiar velocity. 
Taking into account Poisson equation ([3]) we get: 

d 2 (p 

= aR(r)5, a = ^7rGp h o, (13) 

In terms of conformal r and comoving y the hydrodynamic equations can be 
written as: 

£ + « + » W = (14) 
d9 d 

-TT + H6 + — (6v) = -aRS, (15) 
or ay 

We proceed to deal with the Fourier coefficients in the symmetric one- 
dimensional case. 

oo 

5(x,r) = 2j b n (r)exp(ik n x), k n = 2im/L, (16) 



6(x,t)= ^2 a n (r)exp(ik n x), v(x,r) = ^2 exp(ik n x). (17) 

n=— 00 n=— 00 n 
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The reverse transformation is: 



b n ( T ) = L- 1 J dxe lk " x 5(x,r), a n (r) = L' 1 J dx e lk " x 6(x, r) (18) 
o o 

We assume a Q = 0, b = at t = then it is easy to see from (I14|I15I) that 
these equalities are fulfilled for all t > 0. 

The equations for the Fourier coefficients take on the form: 

, oo 

^ + H(r)a n + aR(r)b n + n V = 0, n = ±1, ±2, ... (19) 

GST z — ' V 

p=—oo 1 

^ + a n + nY^^ = 0, n = ±1, ±2,... (20) 
dr ' p 

p=— oo 

The numerical solution of the equations (fl9l) . ( 120|) was performed using the 
4-th order Runge-Kutta method after a truncation of the infinite chain of 
coefficients a n ,b n . Calculations were carried out by the specially written 
GPGPU code using OpenCL SDK by AMD. The time, required to calculate 
a n and b n , in case of 256 values of n (points over k) for single realization of 
initial conditions, is 10 seconds. This is considerably faster than direct using 
of the semi-analytical solution of Section HJ It is important to note that it is 
easy to extend the corresponding algorithms to the 3-D case. 

We calculated coefficients b n of the density contrast as a function of t by 
means of analytical and numerical methods with the same initial conditions. 
On the Fig. [T]these coefficients are presented for t = 0.9 and t = 1.7 for both 
methods. For larger t we observe an infinite growth (for finite time) that 
corresponds to collapse of some plane gravitating shells due violation of con- 
dition (jSJ). Correspondingly, the difference between two methods, that reflects 
the error of calculation, increases for greater t and greater n. For example, in 
order to look how the perturbation propagates from small wavenumbers to 
larger ones, we considered the following initial conditions: b n (0) = where 
integer n varies —128 to 128 except n = ±1; 6±i(0) = 6±i/2 = 0.1; all 
On(0) = 0. For t = 0.2 or t = 0.5 the difference between the values calculated 
by different methods is less than 1%, and for t = 1.7 the difference changes 
from 1.5% to 6.7% as k n = 27m/ L increases from 0.6 to 6. Larger /c-interval 
is presented on Fig. [2J The next figure presents the power spectrum obtained 
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by averaging of the solutions for the ensemble of initial data with uniform 
distribution of 6±i(0), < 6^(0) >= 0.5 . We observe the growth of disper- 
sion, which is explained as follows: as t grows, some of realizations of the 
ensemble of solutions (with larger |6±i(0)|) enter the region which is close to 
the singularity. 

4 Conclusions 

We present a new approach to investigation of the cosmological inhomo- 
geneity by means of the direct integration of hydrodynamic equations in the 
Fourier space. At the moment we studied a one-dimensional hydrodynamical 
evolution of cold (pressureless) gravitating matter. The numerical integration 
has been fulfilled for a random ensemble of initial conditions with subsequent 
averaging procedure to get the power spectrum of the density contrast. We 
used the GPGPU instead of the classic CPU because the problem can be 
easily processed in parallel environment. 

The numerical simulations have been tested using the analytic solution 
that describes the one-dimensional collapse of gravitating shells. The density 
contrast shows a propagation of perturbations from smaller wavenumbers to 
larger ones. The evolution in time ends with a singular growth of the density 
contrast. Correspondingly, we point out a significant growth of dispersion of 
the power spectrum in the non-linear region. 

We consider our results as a first step to the simulations of cosmological 
inhomogeneity growth in the cold matter that could be an alternative to 
the cosmological N-body simulations. We expect that our approach will be 
especially effective in the weakly nonlinear regime. The next step will be 
an implementation of the three-dimensional case of the problem, which is 
technically similar to the one dimensional problem. The trial runs of our 
method allow us to think that it could be really used for power spectrum 
calculations in the 3-D case with realistic requirements to the computer time. 
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Figure 1: Coefficients b n (t) for t = 0.9 (top) and t = 1.7 (bottom) determined 
by semi-analytical and numerical methods with the only nonzero initial values 
6±i(0)=0.1. 
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Figure 2: The evolution of perturbations over k for the initial conditions as 
described on Fig. [TJ 
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Figure 3: Power spectrum <b^> calculated for 100 realizations of randomly 
generated independent initial conditions. 
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